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ABSTRACT 

We present SIGAME simulations of the [Cll] 157.7 /im fine structure line emission from cosmological 
smoothed particle hydrodynamics (SPH) simulations of seven main sequence galaxies at ^ = 2. Using 
sub-grid physics prescriptions the gas in our simulations is modeled as a multi-phased interstellar 
medium (ISM) comprised of molecular gas residing in giant molecular clouds, an atomic gas phase 
associated with photo-dissociation regions (PDRs) at the cloud surfaces, and a diffuse, ionized gas 
phase. Adopting logotropic cloud density profiles and accounting for heating by the local FUV radi¬ 
ation field and cosmic rays by scaling both with local star formation rate (SFR) volume density, we 
calculate the [Cll] emission using a photon escape probability formalism. The [Cll] emission peaks in 
the central < 1 kpc of our galaxies as do the SFR radial profiles, with most [Cll] (> 70%) originating in 
the molecular gas phase, whereas further out (> 2 kpc), the atomic/PDR gas dominates (> 90%) the 
[Cll] emission, no longer tracing on-going star formation. Throughout, the ionized gas contribution 
is negligible (^ 3%). The [Cll] luminosity vs. SFR ([Cll]-SFR) relationship, integrated as well as 
spatially resolved (on scales of 1 kpc), delineated by our simulated galaxies is in good agreement with 
the corresponding relations observed locally and at high redshifts. In our simulations, the molecular 
gas dominates the [Cll] budget at SFR > 20 Mq yr“^ (^sfr ^ 0.5 Mq yr“^ kpc“^), while atomic/PDR 
gas takes over at lower SFRs, suggesting a picture in which [Cll] predominantly traces the molecu¬ 
lar gas in high-density/pressure regions where star formation is on-going, and otherwise reveals the 
atomic/PDR gas phase. 

Subject headings: galaxies: high-redshift - galaxies: ISM - galaxies: star formation - ISM: lines and 
bands 


1. INTRODUCTION 

Single ionized carbon (Cll) can be found throughout 
the interstellar medium (ISM) of galaxies where gas is ex¬ 
posed to UV radiation with energies above the ionization 
potential of neutral carbon (11.3 eV cf. 13.6 eV for hydro¬ 
gen). Cll is found both in regions of ionized and neutral 
gas where, depending on the gas phase, its fine structure 
line [Cii] ^Ps /2 P 1/2 (Vest = 157.714 yum) is collision- 

ally excited by electrons. Hi or H 2 . The ^P 3/2 upper level 
lies 91 K (= hvjk^) above the ^Pi /2 ground state and, 
over a large temperature range (^ 20 —8000 K), the crit¬ 
ical density of [Cll] is only ^ 5 — 44, ^ 1600 — 3800 and 
^ 3300 — 7600 cm“^ for collisions with e“. Hi and H 2 , 
respectively (Goldsmith et al. 2012). As a consequence 
of these characteristics, [Cll] is observed to be one of the 
strongest cooling lines of the ISM, with a line luminosity 
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equivalent to ^ 0.1 — 1% of the far-infrared (FIR) lumi¬ 
nosity of galaxies (e.g., Stacey et al. 1991; Brauher et al. 
2008; Casey et al. 2014). 

Due to high atmospheric opacity at FIR wavelengths, 
observations of [Cll] in the local Universe must be done 
at high altitudes or in space. Indeed, the very first de¬ 
tections of [Cii] towards Galactic objects (Russell et al. 
1980; Stacey et al. 1983; Kurtz et al. 1983) and other 
galaxies (Crawford et al. 1985; Stacey et al. 1991; Mad¬ 
den et al. 1992) were done with airborne observatories 
such as the NASA Lear Jet and the Kuiper Airborne 
Observatory. The Infrared Space Observatory (ISO) al¬ 
lowed for the first systematic [Cll] surveys of local galax¬ 
ies (e.g., Malhotra et al. 1997; Luhman et al. 1998, 2003). 
Detections of [Cll] at high redshifts {z > 1) have also 
become feasible in recent years, with ground-based facil¬ 
ities (e.g., Maiolino et al. 2005, 2009; Hailey-Dunsheath 
et al. 2010; Stacey et al. 2010) and the Herschel Space 
Observatory (e.g., Gullberg et al. 2015). The Atacama 
Large Millimeter Array (ALMA), owing to its tremen¬ 
dous collecting area and high angular resolution, is now 
resolving [Cll] in high-z galaxies (De Breuck et al. 2014; 
Wang et al. 2013) and pushing [Cll] observations of high- 
2 : galaxies to much lower luminosity than before (Ouchi 
et al. 2013; Maiolino et al. 2015; Capak et al. 2015). 

In spite of the observational successes, the interpre¬ 
tation of the [Cii] line as a diagnostic of the ISM and 
the star formation conditions in galaxies is complicated 
by the fact that the [Cll] emission can originate from 
different phases of the ISM. In our Galaxy, about 30% 
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of the total [Cll] emission is found to come from dense 
photo-dominated regions (PDRs), 25 % from cold Hi gas, 
25% from CO-dark H 2 gas, and 20% from ionized gas 
(Pineda et al. 2014). We expect these percentages to 
be different in other galaxies where high levels of star 
formation and/or accretion onto the supermassive black 
hole can boost the energy injected into the ISM and 
change the carbon ionization balance. Models predict 
a reduction in the [Cll] emission from the extreme X-ray 
dominated regions (XDRs) associated with active galac¬ 
tic nuclei (AGNs) (Meijerink et al. 2007; Herrera-Camus 
et al. 2015). In regions of intense FUV-fields (^ 10 x 
the local Galactic field), the CO emitting surfaces of 
clouds will shrink due to photo-dissociation of CO, re¬ 
sulting in a thicker envelope layer of (self-shielding) H 2 
where carbon exist only as Cl or Cll (e.g., Wolfire et al. 
2010). Also, extreme cosmic ray (CR) energy densities 
(i.e., ^ 1000 X that of the Galactic average, as is found in 
local starbursts galaxies; Acero et al. 2009), might lead 
to the efficient destruction of CO (in reactions with He+ 
which forms via cosmic ray ionization), and leave behind 
Cll deep in the FUV-shielded regions of the molecular 
medium (e.g., Bisbas et al. 2015). 

The sensitivity of [Cll] to the presence of FUV radia¬ 
tion led to the line being suggested as a tracer of recent 
star formation. Observations seem to bear out this no¬ 
tion, with normal star-forming galaxies (Lir < 10^^ Lq) 
in the local universe exhibiting a nearly linear relation 
between [Cll] luminosity [L^cii]) arid IR luminosity (Air) 
over several dex (e.g., Malhotra et al. 1997, 2001). Simi¬ 
lar [Cll] — SFR calibrations for local star-forming galax¬ 
ies based on alternative star formation rate (SFR) trac¬ 
ers such as Ha and the UV have also been established 
(e.g., Boselli et al. 2002; de Looze et al. 2011). Recent 
Herschel studies of local spirals and dwarf galaxies have 
found that [Cll] remains a robust tracer of star forma¬ 
tion over scales ranging from ^ 20pc to ^ 1 kpc (e.g., 
De Looze et al. 2014; Herrera-Camus et al. 2015; Ka- 
pala et al. 2015, hereafter referred to as D14, H15, and 
K15, respectively). Pineda et al. (2014) found that our 
Galaxy matches these resolved extragalactic [Cll] — SFR 
relations, both in terms of the slope and overall normal¬ 
ization, only when the [Cll] emission of all the gas phases 
in our Galaxy are combined. 

The scatter in the [Cll] — SFR relation is non-negligible 
- about 0.2 — 0.3 dex in the surface density relations and 
^ 0.3 dex in the luminosity relations (at z ^ 0.3) - and 
must be characterized as best as possible in order to op¬ 
timize [Cll] as a star formation indicator. While some 
of the scatter can be ascribed to imperfect star forma¬ 
tion tracers, it does tend to increase systematically in 
regions of low metallicity, warm dust temperatures, and 
large filling factors of ionized, diffuse gas (D14, H15). 
Understanding this behavior is of importance for studies 
of the low-metallicity, UV-intense environments that we 
expect to encounter in normal galaxies at 2 ; ^ 6, where 
in some cases the [Cll] line is strongly detected (Capak 
et al. 2015), in accordance with the notion that it should 
be one of the brightest available diagnostics of the ISM 
at those epochs, while in other instances the absence of 
[Cll] down to faint flux levels places the sources ^ 10 
times below the local [Cll] — SFR relation (Walter et al. 
2012; Kanekar et al. 2013; Ouchi et al. 2013; Schaerer 


et al. 2015; Maiolino et al. 2015). 

With the [Cll] line increasingly being used as a tracer 
of gas and star formation at high redshifts, efforts have 
also recently been made to simulate the [Cll] emission 
from galaxies. Various sub-grid ‘gas physics’ approaches 
have been applied to both semi-analytical (Popping et al. 
2014a; Munoz & Furlanetto 2014) and hydro dynamical 
(Nagamine et al. 2006; Vallini et al. 2013, 2015) simu¬ 
lations of galaxies. The simulations by Nagamine et al. 
(2006) focus on the [Cll] detectability of lyman break 
galaxies (LBGs) at 2 ; = 3, while Vallini et al. (2013) ap¬ 
ply their simulations to observed upper limits on the [Cll] 
emission from the ^ = 6.6 Lya emitter (LAE) Himiko 
(Ouchi et al. 2013) finding that its metallicity must be 
subsolar. Both set of simulations consider a two-phase 
ISM consisting of a cold and a warm neutral in pres¬ 
sure equilibrium, and both find that the [Cll] emission 
is dominated by the cold neutral medium (CNM). In an 
update of their 2013 simulation, that includes the [Cll] 
contribution from PDRs and a non-uniform metallicity 
distribution in the gas phase, Vallini et al. (2015) finds 
that most of the [Cll] emission from their 2 ; = 6.6 models 
originates in the PDRs with only ^ 10% coming from 
the CNM. 

In this paper we present an adapted version of our code 
Simulator of GAlaxy Millimeter/submillimeter Emission 
(SIGAME; Olsen et al. 2015, submitted) that is capable of 
incorporating [Cll] emission into smoothed particle hy¬ 
drodynamics (SPH) simulations of galaxies. We consider 
a multi-phased ISM consisting of molecular clouds, whose 
surface layers are stratified by EUV-radiation from local¬ 
ized star formation, embedded within a neutral medium 
of atomic gas. In addition, we include the diffuse ionized 
gas intrinsic to the SPH simulations as a third ISM phase. 
The temperatures of the molecular and atomic gas are 
calculated from thermal balance equations sensitive to 
the local EUV-radiation and CR ionization rate. We 
apply SIGAME to GADGET-3 cosmological SPH simula¬ 
tions of seven star-forming galaxies on the main-sequence 
(MS) at 2 ; = 2 in order to simulate the [Cll] emission from 
normal star-forming galaxies at this epoch, examine the 
relative contributions to the emission from the molecular, 
atomic and ionized ISM phases, and the relationship to 
the star formation activity in the galaxies. Throughout, 
we adopt a flat cosmology with Um = 0.27, Ua = 0.73, 
and h = 0.71 (Spergel et al. 2003). 

2. METHODOLOGY OVERVIEW 

SIGAME is applied at the post-processing stage of a SPH 
simulation and takes as its input the following quan¬ 
tities associated with each SPH particle: the position 
{[x,y,z]), velocity {[v^^Vy^Vz]), smoothing length (h), 
gas mass (mspn), hydrogen density (nn), gas kinetic 
temperature (Tk), electron fraction \xq = ng/nn), SER, 
metallicity (Z) as well as the relative abundances of car¬ 
bon ([C/H]) and oxygen ([0/H]). The key steps involved 
in the post-processing are illustrated in Eig. 1 and briefly 
listed below (with details given in subsequent sections): 

1. The SPH gas is separated into its neutral and ion¬ 
ized constituents as dictated by the electron frac¬ 
tion provided by the GAD GET-3 simulations. 

2. The neutral gas is divided into giant molecular 
clouds (GMCs) according to the observed mass 
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Fig. 1. — Schematic illustrating the sub-grid procedures applied to the SPH simulation in post-processing. Each SPH particle is a hybrid 
of neutral and ionized gas. The neutral gas associated with each SPH gas particle is divided into GMCs with masses and sizes following 
the Galactic mass-spectrum and mass-size relation for GMGs. Each GMG has an onion-layer structure, set by the stratification of the 
impinging FUV-field, which consists of outer layer of PDR/atomic gas of Hi and Gll, and an inner molecular region where carbon is found 
in its single ionized state, and in neutral form further in (Section 4.1). The ionized gas associated with each SPH particle is assumed to 
reside in spherical clouds with radii and temperatures given by the SPH smoothing length and gas temperature (Section 4.2). 


function of GMCs in the Milky Way (MW) and 
nearby quiescent galaxies. The GMCs are modeled 
as logotropic spheres with their sizes and internal 
velocity dispersions derived according to pressure- 
normalized scaling relations. 

3. Each GMG is assumed to consist of three spheri¬ 
cally symmetric regions: (1) a FUV-shielded molec¬ 
ular core region where all carbon is locked up in 
Cl and CO; (2) an outer molecular region where 
both Cl and Cll can exist; and (3) a largely neu¬ 
tral atomic layer of Hi, Hll and Cll. The last region 
mimics the FUV-stratified PDRs observed at the 
surfaces of molecular clouds (Hollenbach & Tielens 
1999). This layer can contain both atomic and ion¬ 
ized gas, but we shall refer to it simply as the PDR 
gas. The relative extent of these regions within 
each cloud, and thus the densities at which they 
occur, ultimately depends on the strength of the 
impinging FUV-field and CR ionization rate. The 
latter are set to scale with the local SFR volume 
density and, by requiring thermal balance with the 
cooling from line emission (from Cll, Ol, and H 2 ), 
determine the temperatures of the molecular and 
atomic gas phases. 

4. The remaining ionized gas of the SPH simulation 
is divided into Hll clouds of radius equal to the 
smoothing lengths, temperature equal to that of 
the SPH simulation and constant density. 

5. The [Cll] emission from the molecular, PDR, and 


diffuse ionized gas is calculated separately and 
summed to arrive at the total [Cll] emission from 
the galaxy. In doing so it is assumed that there 
is no radiative coupling between the clouds in the 
galaxy. 

The SPH simulations used in this paper, and the galax¬ 
ies extracted from them, are described in the following 
section. 

3. SPH SIMULATIONS 

Our simulations are evolved with an updated ver¬ 
sion of the public GADGET-3 cosmological SPH code 
(Springel 2005 and S. Huang et al. 2015 in prepara¬ 
tion). It includes cooling processes using the primordial 
abundances as described in Katz et al. (1996), with addi¬ 
tional cooling from metal lines assuming photo-ionization 
equilibrium from Wiersma et al. (2009). We use the 
more recent ‘pressure-entropy’ formulation of SPH which 
resolves mixing issues when compared with standard 
‘density-entropy’ SPH algorithms (see Saitoh & Makino 
2013; Hopkins 2013, for further details). Our code ad¬ 
ditionally implements the time-step limiter of Saitoh & 
Makino (2009), Durier & Dalla Vecchia (2012) which im¬ 
proves the accuracy of the time integration scheme in 
situations where there are sudden changes to a parti¬ 
cle’s internal energy. To prevent artificial fragmentation 
(Schaye & Dalla Vecchia 2008; Robertson & Kravtsov 
2008), we prohibit gas particles from cooling below their 
effective Jeans temperature which ensures that we are al¬ 
ways resolving at least one Jeans mass within a particle’s 
smoothing length. This is very similar to adding pressure 
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TABLE 1 

Global properties of the seven simulated galaxies used for this work at z = 2. 



Gl 

G2 

G3 

G4 

G5 

G6 

G7 

M* [10“ Mq ] 

0.36 

0.78 

0.95 

1.80 

4.03 

5.52 

6.57 

Mg^ [IQio Mq ] 

0.42 

0.68 

1.26 

1.43 

2.59 

2.16 

1.75 

Mneutral [lO^O Mq ] 

0.09 

0.13 

0.57 

0.20 

0.29 

0.29 

0.39 

Mionized [10l° Mq ] 

0.33 

0.55 

0.69 

1.23 

1.30 

1.87 

1.36 

SFR [Mq yr-1 ] 

4.9 

10.0 

8.8 

25.1 

19.9 

59.0 

37.5 

SsFR [Mq yr 1 kpc 

0.016 

0.032 

0.028 

0.080 

0.063 

0.188 

0.119 

Z' 

0.43 

0.85 

0.64 

1.00 

1.00 

1.67 

1.72 


All quantities have been calculated at 2 : = 2 using a fixed cut-out radius of Rent = 10 kpc, which is the radius at which the accumulative 
stellar mass function of each galaxy fiattens. Mgas is the total gas mass, and Mneutral cind Mionized Ihe gas masses in neutral and ionized 
form, respectively (see Section 4). The metallicity (Z' = ZIZq') is the mean of all SPH gas particles within Rent- 


to the ISM as in Springel & Hernquist (2003), Schaye & 
Dalla Vecchia (2008), except instead of directly pressur¬ 
izing the gas we prevent it from cooling and fragmenting 
below the Jeans scale. 

We stochastically form stars within the simulation 
from molecular gas following a Schmidt (1959) law with 
an efficiency of 1% per local free-fall time (Krumholz 
& Tan 2007; Lada et al. 2010). The molecular con¬ 
tent of each gas particle is calculated via the equilib¬ 
rium analytic model of Krumholz et al. (2008, 2009); 
McKee & Krumholz (2010). This model allows us to 
regulate star formation by the local abundance of H 2 
rather than the total gas density, which confines star 
formation to the densest peaks of the ISM. Further im¬ 
plementation details can be found in Thompson et al. 
(2014). Galactic outflows are implemented using the hy¬ 
brid energy/momentum-driven wind (ezw) model fully 
described in Dave et al. (2013); Ford et al. (2015). We 
also account for metal enrichment from Type II super¬ 
novae (SNe), Type la SNe, and AGB stars as described 
in Oppenheimer & Dave (2008). 

3.1. SPH simulations of z = 2 MS galaxies 

We use the cosmological zoom-in simulations presented 
in Thompson et al. (2015), and briefly summarized here. 
Initial conditions were generated using the MUSIC code 
(Hahn & Abel 2011) assuming cosmological parameters 
consistent with constraints from the Planck (Planck Col¬ 
laboration et al. 2014) results, namely = 0.3, Qa = 
0.7, iJo = 70,(78 = 0.8, Us = 0.96. Six target halos 
were selected at z = 2 from a low-resolution Wbody 
simulation consisting of 256^ dark-matter particles in a 
(16h“^Mpc)^ volume with an effective co-moving spatial 
resolution of e = 1.25 h~^ kpc. Each target halo is popu¬ 
lated with higher resolution particles at z = 249, with the 
size of each high resolution region chosen to be 2.5 times 
the maximum radius of the original low-resolution halo. 
The majority of halos in our sample are initialized with a 
single additional level of refinement (e = 0.625 h“^kpc), 
while the two smallest halos are initialized with two ad¬ 
ditional levels of refinement (e = 0.3125 h“^kpc). 

The six halos produce seven star-forming galaxies at 
z = 2 that are free from all low-resolution particles within 
the virial radius of their parent halo. Their stellar masses 
{M^) range from 3.6x 10^ to 6.6x 10^^ Mq and their SFRs 
from 5 to 60 Mq yr“^ (Table 1). We hereafter label the 
galaxies Gl, ..., G7 in order of increasing M^. Other 
relevant global properties directly inferred from the SPH 
simulations, such as total gas mass (Mgas), neutral and 
ionized gas masses (Mneutral and Mionized, respectively). 



Fig. 2. — The SFR-M^c relation at 2; ~ 2 as determined by Spea- 
gle et al. (2014) (dashed line) with the location of our seven simu¬ 
lated galaxies highlighted (filled circles). Dotted-dashed and dot¬ 
ted lines indicate the Icr and 3cr scatter around the relation of Spea- 
gle et al. (2014). For comparison we also show the locus defined by 
3754 1.4 < 2: < 2.5 galaxies from the NEWFIRM Medium-Band 
Survey (gray filled contours), with masses and SFRs calculated 
using a Kroupa IMF (Whitaker et al. 2011). 

average SFR surface density (Esfr), and average metal- 
licity {Z')^ can also be found in Table 1. 

Fig. 2 shows the locations of Gl, ..., G7 in the SFR—M^ 
diagram. The galaxies are consistent with observational 
determinations of the z ^ 2 MS of star-forming galaxies 
(Whitaker et al. 2011; Speagle et al. 2014). 

4. MODELING THE ISM 

As illustrated in Fig. 1, the first step in modeling the 
ISM is to split each SPH particle into an ionized and a 
neutral gas component. This is done using the electron 
fraction, Xg, associated with each SPH particle, i.e.,: 

^neutral (I (I) 

^ionized ~ ^e^SPYL- (2) 

The electron fraction from GADGET-3 gives the den¬ 
sity of electrons relative to that of hydrogen, ng/nn, and 
can therefore reach values of ^ 1.16 in the case where 
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helium is also ionized. As a result we re-normalized the 
distribution of values to a maximal of 1 so as to 
not exceed the total gas mass in the simulation. Fig. 
3 shows the distribution of SPH gas particles masses in 
G4 - chosen for its position near the center of the stel¬ 
lar and gas mass ranges of Gl, ..., G7 - along with the 
mass distributions of the neutral and ionized gas com¬ 
ponents obtained from eqs. 1 and 2. The ionized gas is 
seen to have a relatively flat distribution spanning the 
mass range ^ 10^-^ — 10^-^ Mq. The neutral gas, how¬ 
ever, peaks at two characteristic masses ( r\j 10^-^ Mq and 
^ 10^-^ Mq), where the lower mass peak represents gas 
particles left over from the first generation of stars in the 
simulation. 



Fig. 3.— The distribution of SPH gas particle masses (solid his¬ 
togram) in G4. The distribution peaks at two characteristic masses 
10^'^ Mq and 10^'® Mq), where the lower mass peak repre¬ 
sents gas particles left over from the first generation of stars in the 
simulation. Splitting the gas into its neutral and ionized compo¬ 
nents according to eqs. 1 and 2 results in the mass distribution 
given by the dotted and dashed histograms, respectively. 


4.1. GMCs 

4.1.1. Masses, sizes and veloeity dispersion 

The neutral gas mass, mneutrab associated with a given 
SPH particle is divided into GMCs by randomly sam¬ 
pling the GMC mass spectrum as observed in the Galac¬ 
tic disk and Local Group galaxies: oc with 

(3 = l.S (Blitz et al. 2007). Similar to Narayanan et al. 
(2008a,b), a lower and upper cut in mass of 10^ Mq and 
10^ Mq, respectively, are enforced in order to ensure the 
GMC masses stay within the range observed by Blitz 
et al. (2007). Up to 40 GMCs are created per SPH par¬ 
ticle but typically most (> 90%) of the SPH particles 
are split into four GMCs or less. While mneutrai never 
exceeds the upper GMC mass limit (10^ Mq), there are 
instances where rUneutrai is below the lower GMC mass 
limit (10^ Mq). In those cases we simply discard the gas, 
i.e., remove it from any further sub-grid processing. For 
the highest resolution simulations in our sample (Gl and 
G2), the discarded neutral gas amounts to ^ 6% of the 
total neutral gas mass, and < 0.02 % in the remaining 
five galaxies. We shall therefore assume that it does not 
affect our results significantly. 


The GMCs are randomly distributed within 0.2 x the 
smoothing length of the original SPH particle, but with 
the radial displacement scaling inversely with GMC mass 
in order to retain the original gas mass distribution as 
closely as possible. To preserve the overall gas kinematics 
as best possible, all GMCs associated with a given SPH 
particle are given the same velocity as that of the SPH 
particle. 

GMC sizes are obtained from the pressure-normalized 
scaling relations for virialized molecular clouds which re¬ 
late cloud radius (i^cMc) with mass (mcMc) and exter¬ 
nal pressure (Pext)- 

Pgmc _ f Pext/fa A f rriGuc 

pc “Vio'^cm-^Ky V^QOMq; ‘ 

We assume Pext = Ptot/(l + cto + /^o) foi* relative cosmic 
and magnetic pressure contributions of cto = 0.4 and 
Po = 0.25 (Elmegreen 1989). For the total pressure (Ptot) 
we adopt the external hydrostatic pressure at mid-plane 
for a rotating disk of gas and stars, i.e.,: 


Ptot 



gas 




(4) 


where Egas and are the local surface densities of gas 
and stars, respectively, and crgas,± and their lo¬ 

cal velocity dispersions measured perpendicular to the 
mid-plane (see e.g., Elmegreen (1989); Swinbank et al. 
(2011)). Eor each SPH particle, all of these quantities 
are calculated directly from the simulation output (using 
a radius of P = 1 kpc from each SPH particle), and it is 
assumed that the resulting Pext is the external pressure 
experienced by all of the GMCs generated by the SPH 
particle. We find that GMCs in our simulated galax¬ 
ies are subjected to a wide range of external pressures 
(Pext/ r\j 10^ — 10^ cm ^ K). Eor comparison, the range 
of pressures experienced by clouds in our Galaxy and in 
Local Group galaxies is Pext/te ^ 10^ —10^ cm“^ K with 
an average of Pext/^B ^ 10^ cm“^ K in Galactic clouds 
(Elmegreen 1989; Blitz et al. 2007). This results in GMC 
sizes in our simulations ranging from Pgmc = 1 — 300 pc. 

The internal velocity dispersion (cr^) of the GMCs is 
inferred from the virial theorem, which provides ns with 
a pressure-normalized ay — Pgmc relation: 


= 1.2 km 


.-1 


/^B A f Pgmc \ 


10^ cm ^ K 


V pc / 


(5) 


where the normalization of 1.2 km s“^ comes from studies 
of Galactic GMCs (Larson 1981; Elmegreen 1989; Swin¬ 
bank et al. 2011). 


4.1.2. GMC density and temperature strueture 

We assume a truncated logotropic profile for the total 
hydrogen number density of the GMCs, i.e.,: 

^h(P) = ^H,ext ^ ^ , (6) 


where nH(P > Pgmc) = 0. Eor such a density profile it 
can be shown that the external density, nH,ext 5 is 2/3 of 
the average density: 


^H,ext = 2/3(nH) = 2/3 


^GMC 

I/Stttoh-Rgmc^ 


(7) 
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While the total hydrogen density follows a logotropic pro¬ 
file, the transition from H 2 —t^Hi/Hii is assumed to be 
sharp. Similarly for the transition from Cl— t^Cii. This 
is illustrated in Fig. 4, which shows an example density 
profile of a CMC from our simulations. From the center 
of the CMC and out to hydrogen is in molecular 

form. Beyond Ru 2 ^ hydrogen is found as Hi and Hll out 
to i^GMC- 



0 5 10 15 20 25 

R [pc] 


Fig. 4.— Example H number density profile (solid line) for a 
GMC of mass mcMC = 1-3 x 10^ Mq and radius -Rqmc = 21 pc. 
Also shown are the density profiles of H 2 (= O.Stih), Hi and Hii. 
Note, the transition from molecular to atomic H is assumed to 
happen instantaneously at Rii 2 • Ttie total C abundance follows 
that of H but scaled with the [C/H] abundance (as provided by 
the ‘ps-rent’ SPH particle, see Section 2). Cll can exist throughout 
the cloud except for the very inner region (R < i^cil indicated in 
brown), with the Cl to Cll transition happening instantaneously 
at Rci- [Cll] emission from the layer Rci < R < R 112 (indicated 
in red) is referred to as ‘molecular’ emission, while [Cll] emission 
from Rii^ < R < Rqmc (indicated in orange) is referred to as 
‘PDR’ emission. The relative thickness of these layers is set by 
the impinging FUV radiation field and cosmic rays (illustrated as 
purple and blue arrows, respectively), with the former undergoing 
attenuation further into the cloud. 

The size of the molecular region, when adopting the 
logotropic density profile, is related to the total molecular 
gas mass fraction (/^oj) of each GMC: 

rriGMC 

!t^RdR y .g. 

\Rgmc) 

’’Ha = \J (10) 

where trs is the fractional radius, = Rn^/R gmc- We 
find by assuming Hi H 2 equilibrium and using 
the analytical steady-state approach of Pelupessy et al. 
(2006) for inferring for a logotropic cloud subjected 
to a radially incident FUV radiation field (see also Olsen 
et al. 2015). In this framework, the value of (and 
thereby trs) depends on: 


1. The cloud boundary pressure (Pext), which is cal¬ 
culated as explained in Section 4.1.1. 

2. The metallicity {Z') oi the GMC, which is inherited 
from its parent SPH particle and assumed constant 
throughout the cloud. 

3. The kinetic temperature of the gas at the GMC 
surface (T/c(i?GMc))- This quantity is calculated in 
an iterative process together with by solving 
the following thermal balance equation: 

FpE + rGR,Hi = Agii + Aqi, (11) 

where Fre is the heating rate associated with the 
photoelectric ejection of electrons from dust grains 
by the FUV field and Fgr^i is the heating rate by 
cosmic rays in atomic gas. The main cooling agents 
are assumed to be due to [Cll] and [Ol] (63 /im and 
145/im) line emission (i.e., Agh and Aqi, respec¬ 
tively). FpE, rGR,Hi, and Agh all depend on the 
electron fraction at i^GMG, which is determined by 
the degree of Hi ionization by the local FUV ra¬ 
diation field and CR ionization rate (see below for 
how these quantities are derived). For analytical 
expressions for the heating rates, we refer to Olsen 
et al. (2015). For Aqi we use the expressions given 
in Rollig et al. (2006). The calculation of Agh at 
the GMC surface is detailed in Section 5 and Ap¬ 
pendix A. 

4. The strength of the local FUV radiation field (Go) 
and the CR ionization rate (Cgr) impinging on the 
GMCs. These quantities do not come out from 
the simulation directly, and instead they are calcu¬ 
lated by scaling the Galactic FUV field (Go,mw) 
and CR ionization rate (Cgr,mw) with the lo¬ 
cal SFR volume density in the simulations, i.e.. 
Go (X Go,MW (SFRDiocai/SFRDMw) and Cgr oc 
CcR,MW (SFRDiocai/SFRDMw), where SFRDiocai 
is estimated for each SPH particle as the vol¬ 
ume averaged SFR within a 2kpc radius. We 
have adopted Milky Way values of Go,mw = 
0.6Habing (Seon et al. 2011) and Ccr,mw = 3 x 
10“^^ s“^ (Webber 1998). For SFRDmw we adopt 
0.0024 Mq yr“^ kpc“^, inferred from the average 
Galactic SFR (O.3M0yr“^) within a disk lOkpc 
in radius and 0.2 kpc in height (Heiderman et al. 
2010; Bovy et al. 2012). 

5. The electron fraction at the cloud boundary. This 
fraction is not the previously introduced Xg, which 
was inherent to the SPH simulations and used to 
split the SPH gas into a neutral and ionized gas 
phase. Instead, it is the electron fraction given by 
the degree of ionization of Hi caused by the Go and 
Cgr impinging on the cloud. This fraction (and 
thus the Hl:Hll ratio) is calculated with CLOUDY 
vl3.03 (Ferland et al. 2013) given the hydrogen 
density and temperature at the cloud boundary 
and assuming an unattenuated Go and Ccr at 
Rgmg- 

The [Cll] emitting region in each of our GMCs is de¬ 
fined as the layer between the surface of the cloud and the 
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depth at which the abundances of C and are equal. 
Hence, if the latter occurs at a radius Rqi from the cloud 
center, the thickness of the layer is Rqmc — Rci (Fig. 
4). At radii < carbon atoms are for simplicity 

assumed to be in neutral form. In order to determine 
the fractional radius rci (= Rci/R gmc) we follow the 
work of Rollig et al. (2006) (but see also Pelupessy & 
Papadopoulos 2009), who considers the following domi¬ 
nant reaction channels for the formation and destruction 
of C+: 

C + 7^C+ + e- (12) 

C+ + e-^C + 7 (13) 

C+ + H2 ^ CH+ + 7. (14) 

In this case, rci can be found by solving the following 
equation: 

g—A<'^FUV^v(^Ci) 

5.13 X lO-^V^Go / -^- dn (15) 

Jl ^ 

= nH(rci)[ac^c + O.bkc], 

where the left-hand side is the C'*" formation rate due 
to photo-ionization by the attenuated FUV field at rci 
(eq. 12), and the right-hand side is the destruction rate 
of C+ due to recombination and radiative association 
(eqs. 13 and 14). The constants ac = 3 x 10“^^ cm“^ 
s“^ and kc = S X 10“^^ cm“^ s“^ are the recombination 
and radiative association rate coefficients. Note, we have 
accounted for an isotropic FUV field since /i = cos6>, 
where 0 is the angle between the Poynting vector and 
the normal direction. Av(rci) is the visual extinction 
corresponding to the Rqmc — Rci layer, and is given 
by Av(rci) = 0.724crdust^'(^H)^GMC In (rcU^), where 
c^dust = 4.9 X 10“^^ cm^ is the FUV dust absorption cross 
section (Pelupessy & Papadopoulos 2009; Mezger et al. 
1982). <^FUV accounts for the difference in opacity be¬ 
tween visual and FUV light and is set to 3.02. Xq is 
the carbon abundance relative to H and is calculated by 
adopting the carbon mass fractions of the parent SPH 
particle (self-consistently calculated as part of the overall 
SPH simulation) and assuming it to be constant through¬ 
out the GMC. 

The temperature of the [Gll]-emitting molecular gas 
(i.e., from the gas layer between Rqi and Rhs) is assumed 
to be constant and equal to the temperature at Ru^ . The 
latter is given by the thermal balance: 

FpE + rcR,H2 = Ah 2 + Acii + Aoi, (16) 

where rcR,H 2 i^ i^Fe CR heating in molecular gas, and 
Ah 2 is the cooling due to the S(0) and S(l) rotational 
lines of H 2 (Papadopoulos et al. (2014); see also Olsen 
et al. 2015). rcR,H 2 depends on Xg, which is calcu¬ 
lated with CLOUDY for the local CR ionization rate and 
the attenuated FUV field at i.e. Goe“'^^uvAv(rH 2 ), 
Av(^H 2 ) corresponds the extinction through the outer 
atomic transition layer of each GMC. 

The temperature of the PDR gas (i.e., from the gas 
between Ru^ and Rgmc) is assumed to be constant and 
equal to the temperature at Rgmg- 

For each GMC we solve in an iterative manner simul¬ 
taneously for Rhs (eq- 10) and Ra (eq. 15), as well as 



iog( >^905 [M©] ) 




Fig. 5. — Top: the number distribution of GMC masses (solid 
histogram) together with the number distributions of the molecu¬ 
lar (dashed histogram) and PDR (dotted histogram) gas compo¬ 
nent. Middle: GMC-mass-weighted distributions of Rci (dotted), 
Rh 2 (dashed) and Rqmc for GMCs in G4. Bottom: GMC-mass- 
weighted distributions of the temperature at (dashed) and at 
the cloud surfaces (solid) for GMCs in G4. 
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for the gas temperature at Rqmc and at i?H 2 (eqs. 11 
and 16, respectively) 

The resulting distributions of and Rci for the 
CMC population in G4 are shown in Fig. 5 (middle 
panel), along with the distribution of Rqmc (obtained 
from eq. 3). In GMCs in general, we expect Rqi < Rb. 2 ^ 
due to efficient H 2 selfshielding. However, as Fig. 5 shows 
some of the GMGs in our simulations have very small 
i?H 2 (due to them having virtually zero molecular gas 
fractions), and in those cases Ra can be equal to or even 
exceed Rn ^. The latter implies that the [Gll] emission is 
only coming from the PDR phase, with no contribution 
from the molecular gas. 

The distributions of the kinetic temperature of the gas 
at i ? H2 Rqmc for the GMG population in G4 are 
also shown in Fig. 5 (bottom panel). The temperatures 
at the cloud surfaces range from ^ 8K to ^ 10^-^ K, and 
the temperatures at Rn^ lies between ^ 8 and 1800 K. 

With i?H 2 und Ra determined for each GMG we can 
calculate the gas masses associated with the molecular 
and PDR gas phase, respectively. The resulting mass 
distributions are shown in the top panel of Fig. 5, along 
with the distribution of total GMG masses {ttigmc) us 
determined by the adopted GMG mass spectrum (Sec¬ 
tion 4.1.1). We see that most of the molecular and PDR 
gas masses follow the total GMG mass spectrum. There 
is, however, a fraction of GMGs with extremely small 
molecular gas masses (corresponding to ^ 0). On 
the other hand, there are some GMGs with small PDR 
gas masses, i.e., clouds that are so shielded from FUV 
radiation that they are almost entirely molecular. 

4.2. The ionized gas 

The ionized gas in our simulations (see eq. 2), is as¬ 
sumed to be distributed in spherical clouds of uniform 
densities and with radii (Run) equal to the smoothing 
lengths of the original SPH particles. These ionized re¬ 
gions in our simulations are furthermore assumed to be 
isothermal with the temperatures equal to that of the 
SPH gas. Fig. 6 shows the size (top) and temperature 
(bottom) distribution for the ionized clouds in G4. Gloud 
sizes range from ^ 0.1 to ^ lOkpc, with more than 
60% of the ionized gas mass residing in clouds of size 
< 1000 pc. For comparison, the range of observed sizes 
of Hll clouds in nearby galaxies is 10 — 1000 pc (Oey & 
Glarke 1997; Hodge et al. 1999). The temperatures range 
from ^ 10^ K to ^ 10^ K with the bulk of the ionized gas 
having temperatures ^ 10^“^ K. 

5. THE [Cii] LINE EMISSION 

The [Gll] luminosity of a region of gas is the volume- 
integral of the effective [Cll] cooling rate per volume, 
i.e.,: 

^[Cii] = [ AciidV. (17) 

JAV 

Since we have adopted spherical symmetry in our sub¬ 
grid treatment of the clouds (both neutral and ionized), 
we have: 

rR2 

-^[Cii] ~ ^TT / KqiiR? dR^ (18) 

jRi 




Fig. 6. — Mass-weighted histograms of the size (top) and tem¬ 
perature (bottom) distributions of the ionized clouds in G4. 

where Ri = Ra and R 2 = Rh 2 foi* Ihe molecular phase; 
Ri = 7?h 2 and R 2 = Rgmc foi* the PDR region, and 
= 0 and R 2 = Rhu for the ionized gas. 

The effective [Gll] cooling rate is: 

Acii = A^ijdf^ncuhn, (19) 

where A^\ (=2.3x10“^ s“^)^ is the Einstein coefficient 
for spontaneous decay. We ignore the effects of any 
background radiation field from dust and the cosmic 
microwave background (GMB). P is the [Gll] photon 
escape probability for a spherical geometry (i.e., (3 = 
(1 — exp{—T))lT^ where r is the [Gll] optical depth), 
is the fraction of singly ionized carbon in the upper ^P ^/2 
level and is determined by radiative processes and colli- 
sional (de)excitation (see Appendix A). The latter can 
occur via collisions with e“. Hi, or H 2 , depending on 
the state of the gas. In our simulations the collisional 
partner is H 2 in the molecular phase. Hi and e~ in the 
PDR regions, and e~ in the ionized gas. Analytical ex¬ 
pressions for the corresponding collision rate coefficients 
as a function of temperature are given in Appendix A. 
ncii is the number density of singly ionized carbon and 
is given by ncn = -Ac/ch^h, where fcn is the fraction of 

^ Einstein coefficient for spontaneous emission taken from the 
LAMDA database: http://www.strw.leidenuniv.nl/~moldata/, 
Schoier et al. (2005) 
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carbon atoms in the singly ionized state. For the latter 
we used tabulated fractions from CLOUDY vl3.03 over a 
wide range in temperature, hydrogen density, FUV field 
strength, and CR ionization rate. 

We calculate the integral in eq. 18 numerically by split¬ 
ting the R 2 — Ri region up into 100 radial bins. In each 
bin, nn is set to be constant and - in the case of the 
molecular and PDR regions - given by the logotropic 
density profile at the radius of the given bin (Fig. 4). For 
the ionized clouds, nn is constant throughout (Section 
4.2). For the PDR regions (i.e., from Rh 2 fo Rgmc) 
we assume that the temperature, electron fraction and 
Go are kept fixed to the outer boundary value at Rgmc 
( i.e., no attenuation of the FUV field). This implies that 
the [Cll] luminosity from the PDR gas is an upper limit. 
Similar for the [Cll] emission from the molecular region 
(i.e., from Rqi to RH 2 ): where we assume that the tem¬ 
perature throughout this region is fixed to its values at 
Rh 2 - Also, throughout this region we adopt the attenu¬ 
ated FUV field at Rh 2 - 

6. RESULTS AND DISCUSSION 

Having divided the ISM in our galaxies into molecu¬ 
lar, atomic and ionized gas phases, and having devised a 
methodology for calculating their [Cll] emission, we are 
now in a position to quantify the relative contributions 
from the aforementioned gas phases to the total [Cll] 
emission, and examine their relationship to the on-going 
star formation. 

6 . 1 . Radial [Cll] luminosity profiles 

First, however, to get a sense of the distribution of gas 
and star formation in our simulated galaxies, we show 
in Fig. 7 surface density maps of the total SPH gas 
(left column) and star formation rate (middle column) 
when viewed face-on. The maps reveal spiral galaxy mor¬ 
phologies, albeit with some variety: some (Gl, G 2 and 
G3) show perturbed spiral arms due to on-going mergers 
with satellite galaxies; others (G4, G5, G 6 and G7) have 
seemingly undisturbed, grand-design spiral arms; a cen¬ 
tral bar-like structure is also seen in some (G 2 , G4, G5 
and G7). Overall, the star formation is seen to be much 
more centrally concentrated than the SPH gas. This is 
especially true for Gl and G 2 , which have very centrally 
peaked star formation. The radial SFR and [Cll] lumi¬ 
nosity profiles of Gl, ..., G7 - derived by summing up 
the SFR and the [Gll] luminosity within concentric rings 
(of fixed width: 0.2 kpc) - are also shown in Fig. 7. We 
have inferred the radial [Cll] luminosity distribution for 
the full ISM as well as for the individual gas phases. 

The radial SFR profiles of our model galaxies typically 
peak at R ^ 0.5 kpc (in Gl at R ^ 0.5 kpc) and then 
tail off with radius. In some cases, local peaks in the 
star formation activity occur at galactocentric distances 
> 2 kpc, corresponding to the locations of either satellite 
galaxies (G 2 and G3) or spiral arms (G5 and G7). 

The total [Gll] luminosity profiles (gray histograms) 
also peak at R ^ 0.5 kpc, and within the central R ^ 
1.5 kpc there is in general a good correspondence between 
the total [Cll] emission and the star formation activity. 
This correspondence is driven by the molecular gas phase 
which dominates the [Cll] emission in the central regions. 
The molecular [Cll] emission is seen to correlate strongly 
with the star formation out to radii ^ 5 kpc and beyond 


(e.g., G3 and G4). There are cases, however, where lo¬ 
calized enhancements in the SFR are not matched by 
increased [Gll] emission from the molecular gas (e.g. G 2 , 
and G3). These SFR enhancements are reflected in the 
[Cll] emission profile of the ionized gas, which at galacto¬ 
centric distances ^ 2 kpc follow the SFR closely (despite 
contributing only a small fraction to the total [Gll] emis¬ 
sion budget, see below). In contrast, the [Cll] emission 
from the PDR gas, which dominates the total [Cll] lu¬ 
minosity at 2.0 kpc < R ^ 10 kpc, does not appear to be 
a sensitive tracer of the SFR. At R ^ 1.5 kpc where the 
SFR peaks, the [Cll] emission from this phase is seen to 
drop. At larger radii, the PDR [Cll] emission declines 
but at a more gradual rate than the star formation. 

In the bottom right panel of Fig. 7 we show the Galac¬ 
tic SFR and [Cll] luminosity radial profiles from Pineda 
et al. (2014), who observed the [Cll] emission from ( 1 ) 
CO-dark H 2 gas and dense PDRs, ( 2 ) cold neutral Hi 
gas, and (3) hot ionized gas in our own Galaxy. In order 
to facilitate an approximate comparison with our simu¬ 
lations we identify these three Galactic ISM phases with 
the molecular, atomic, and ionized gas in our simulations. 
It is important to keep in mind, however, that the ISM 
in our simulations has a higher pressure, is kinematically 
more violent, and is more actively forming stars than is 
the case in our Galaxy. 

The SFR and [Cll] profiles in our Galaxy peak at larger 
radii {R ^ 4 —5 kpc) than in our simulated galaxies. This 
is not surprising given the low content of star formation 
and gas in the Galactic bulge, and the fact that the bulk 
of star formation in our Galaxy takes place in the disk. 
In contrast, gas is still being funneled toward the central 
regions of our simulated galaxies where it is converted 
into stars. Thus the SFR level in our simulated galaxies 
is much higher (by ^ 10 x) and more centrally concen¬ 
trated than in our Galaxy, where stars form at a rate of 
< O.I — IM 0 yr“^ across the disk. Remarkably, signifi¬ 
cant levels of [Cll] emission extend out to R ^ 20 kpc in 
our Galaxy, well beyond the point where star formation 
has ceased. Here, the emission is completely dominated 
by CO-dark H 2 gas and dense PDR regions. This is sim¬ 
ilar to the picture seen in our simulated galaxies where 
the [Cll] emission at large radii is dominated by a neu¬ 
tral gas phase - designated PDR gas in our simulations, 
dubbed CO-dark H 2 + dense PDR gas in Pineda et al. 
(2014) - that is largely uncoupled from star formation. 
The radial [Cll] profile of the ionized gas in the Galaxy 
roughly follows the SFR profile, as is the case in our 
simulations. 

Fig. 8 displays the fractional [Cll] luminosity from the 
different ISM phases: top panel for the entire disk {R < 
10 kpc), middle panel for the central region {R < I kpc), 
and bottom panel for the outer disk {R > 2 kpc). Within 
R < 10 kpc, the molecular gas can constitute from ^ 31% 
(G 2 ) to ^ 91 % (G7) of the total [Cll] luminosity; for the 
PDR gas the range is ^ 9% (G7) to ^ 67% (G 2 ). Fig. 
8 shows that the contribution from the molecular gas to 
the total [Cll] emission increases with the overall SFR 
of the galaxy. A reverse trend is seen for the PDR gas. 
As expected, the total [Cll] emission from the central 
regions {R < I kpc) is dominated by the molecular gas (> 
70%) while further our {R > 2 kpc) the PDR gas phase 
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Fig. 7.— Gas and SFR surface density maps (left and middle columns, respectively) of our SPH simulated galaxies viewed face-on (Gl, 
G7 from top to bottom; see also Thompson et al. (2015)). The horizontal white bars correspond to a physical scale of 5kpc. The 
right-hand column shows the radial profiles of the total [Gll] luminosity (gray histogram) and the contributions from molecular gas (red 
curve), PDR (orange curve) and ionized gas (blue curve). The SFR radial profiles are also shown (black curve). The radial profiles were 
determined by summing up the [Gll] luminosity and SFR within concentric rings with fixed width of 0.2 kpc. 
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Fig. 7. — Continued. In the bottom panel we show for reference the radial profile of the Galaxy (gray histogram) along with the 

contributions from CO-dark H 2 gas and dense PDRs (red curve), cold Hi (orange curve), and ionized gas (blue curve) (Pineda et al. 2014). 
Also shown is the radial SFR profile (black curve) inferred from the 1.4 GHz intensity distribution. Fixed radial bin widths of 0.5 kpc have 
been adopted (see Pineda et al. (2014) for details). 


dominates (> 90%) (see Fig. 8, bottom two panels). For 
reference, we note that in our own Galaxy about 55 % of 
the total Galactic [Gll] luminosity (within R < 20 kpc) is 
from molecular gas and dense PDRs, 25 % from cold Hi, 
and 20 % from the ionized gas (Pineda et al. 2014). Thus, 
the ionized phase is a more important contributor to the 
overall [Gll] budget in our Galaxy than in the simulated 
galaxies presented here where the contribution from the 
ionized gas is ^ 3% in all cases. 

6.2. The integrated L^cu] ~ SFR relation 

Fig. 9 shows the integrated L^cu] ~ SFR relations for 
our simulated galaxies: top panel for the full ISM and, 
in separate panels below, for each of the three ISM phases 


considered in our simulations. 

When considering the entire ISM a tight correlation 
between L^qu] SFR emerges, which is well fit in log- 
log space by a straight line with slope 1.27 ± 0.17 (solid 
line in the top panel). This relation is largely set by 
the molecular and PDR gas phases. The molecular gas 
phase, itself exhibiting a strong correlation between [Cll] 
and SFR with slope 1.72T0.22, drives the slope of the to¬ 
tal correlation. The PDR gas on the other hand, shows a 
weaker [Cll] —SFR correlation with a slope of 0.43 ±0.20 
but contributes significantly to the normalization of the 
total [Cll] — SFR relation, especially at the low SFR end 
(see Fig. 9). The [Gll] emission from the ionized gas also 
shows a weak dependency on SFR (slope 0.44±0.30) but 
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R < lOkpc 



Fig. 8.— Contributions to the total [Cll] luminosity from the 
molecular (red), PDR (yellow), and ionized (blue) gas phases for 
each of our simulated galaxies, ranked according to their total 
SFRs. The three panels show the relative contributions within 
R < lOkpc (top), for R < 1 kpc (middle), and for R > 2kpc 
(bottom). The horizontal dashed lines and the percentages given 
indicate the relative contributions to the total [Cll] luminosity of 
our own Galaxy from CO-dark H 2 + dense PDR gas (red; 55%), 
cold atomic gas (yellow; 25%), and ionized gas (blue; 20%) (Pineda 
et al. 2014). The total SFR of the Galaxy is 1.9 Mq yr“^ (Chomiuk 
& Povich 2011). 

does not contribute significantly to the total [Cll] emis¬ 
sion, its normalization factor being > lOx below that of 
the molecular and PDR gas (Fig. 9, bottom panel). 

In the top panel of Fig. 9 we compare the L^qu] ~ SFR 
relation obtained from our simulated galaxies with sam¬ 
ples of [Cll] detected galaxies in the redshift range ^ ^ 
0—6.5 compiled from the literature. Our simulated galax¬ 
ies are seen to match the observed L^qu] ~ SFR relation 
both in terms of the slope of the relation and its overall 
normalization. A power-law fit to the simulated galaxies 
(shown as solid line in Fig. 9) yields a near-linear slope 
(1.27±0.17). Normal star-forming galaxies at ^ ^ 0, with 
similar levels of SFR (^ 0.2 — lOOMQyr”^; Malhotra 
et al. 2001) as our simulated galaxies, are consistent with 
a linear correlation given by L^cu] = 1-05 x lO^SFR (dot¬ 
ted line in Fig. 9) with a scatter of 0.3 dex (Magdis et al. 
2014)^. The scatter of our simulated galaxies around 
their best-fit relation is 0.14 dex, i.e., significantly lower. 
We attribute this to the fact that our simulated galax¬ 
ies constitute a fairly homogeneous (and small) sample 
spanning a rather small range in SFR, L^cu]^ Z', un¬ 
like the observed samples with which we are comparing. 

A direct comparison with [Cll]-detected galaxies at 
high redshifts is complicated by the fact that the latter 

® This expression is inferred from a power-law fit by Magdis 
et al. (2014) to the [Gll] and FIR (42.5 — 122.5 p-m) luminosities 
(in units of Lq) of the Malhotra et al. (2001) sample: = 

10-2.51±o.392^^j^, have made use of the conversion 

SFR = Lfir/3.4 X 10^. 


typically have significantly larger SFRs than our model 
galaxies. Furthermore, high- 2 ; samples are often hetero¬ 
geneous (e.g., significant AGN contribution, cf. Gullberg 
et al. 2015). One exception, however, is the recent [Cll]- 
detected sample of star-forming galaxies at ^ ^ 5 — 6 pre¬ 
sented by Capak et al. (2015) (shown as magenta stars 
in Fig. 9), which span the same range in SFR as our 
simulations. These galaxies are seen to be in excellent 
agreement with the L^qu] ~ SFR relation defined by our 
simulations, both in terms of slope and normalization. 

Extrapolating our best-fit relation to SFRs > 
300 Mq yr“^ in order to compare with other high- 2 : sam¬ 
ples, the relation is seen to overshoot the data. A 
power-law fit to the z > 1 star-forming galaxies with 
SFR> 300 Mq yr“^ compiled by D14 yields: L^cu] = 

1.7 X lO^SFR^’^^ (D14; shown as the dashed line in Fig. 
9), i.e., formally, a shallower relation than that of our 
simulated galaxies and that of the 2 : ^ 0 sample (al¬ 
beit less so). Finally, we stress that the [Cll]-detected 
galaxies at ^ > 1 with SFRs ^ 300MQyr“^ likely de¬ 
rive from rather complex environments (e.g. Narayanan 
et al. 2015), which may not correspond to the relatively 
quiescent MS star-forming galaxies modeled here. 

6.3. The resolved E[cn] ~ ^sfr relation 

In Fig. 10 we show the combined ~ ^sfr relation 
of all seven simulated galaxies in their face-on configu¬ 
ration. The relation is shown for the entire ISM (top 
panel) and for each of the separate gas phases (bottom 
three panels). Surface densities were determined within 
1 kpc X 1 kpc regions. Contours reflect the number of re¬ 
gions at a given (IlsFR 5 ^[Cii])-combination and are given 
as percentages of the peak number of regions. 

Given the variations in the local star formation condi¬ 
tions within our model galaxies, we can explore the re¬ 
lationship between [Cll] and star formation over a much 
wider range of star formation intensities than is possible 
with the integrated quantities. A T>^cu] ~ ^sfr correla¬ 
tion spanning more than five decades in Esfr is seen for 
the entire ISM as well as for the individual ISM phases. 
Similar to the integrated [Cll] vs. SFR relations in the 
previous section, the molecular and PDR phases domi¬ 
nate the resolved [Cll] emission budget at all SFR surface 
densities, and with the molecular relation being steepest 
and exhibiting the largest degree of scatter. The resolved 
[Cll] emission from the ionized phase, while clearly corre¬ 
lated with EsfR: is largely negligible at all star formation 
densities. 

We compare the — Esfr relation for our simu¬ 

lated galaxies to that of three resolved surveys of nearby 
galaxies: (1) the ^ 50pc-scale relation derived from five 
3' X 3' fields toward M31 (K15; shown as green points 
in Fig. 10), (2) the kpc-scale relation obtained for 48 
local dwarf galaxies covering a wide range in metallic- 
ities {ZjZ^ = 0.02 — 1, D14; cyan crosses), and (3) 
the kpc-scale relations for local, mostly spiral, galaxies 
(H15; magenta contours). The data from these surveys 
have been converted to the Chabrier IMF assumed by 
our simulations by multiplying Esfr with a factor 0.94 
when a Kroupa IMF was adopted (D14; K15) or 0.92 
when a truncated Salpeter IMF was used (H15), fol¬ 
lowing Galzetti et al. (2007) and Speagle et al. (2014). 
Regarding the observations by HI5, we adopt the raw 
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Fig. 9. — vs. SFR for our simulated galaxies (big filled circles). The same color-coding as in Fig. 2 is used and listed again for 

convenience in the bottom panel. From top to bottom panel we show the [Cll] luminosity from the full ISM (gray), the molecular gas (red), 
PDRs (orange), and ionized gas (blue). For comparison, we show individual [Cll] observations of 54 2 ; ~ 0 spirals (purple filled circles; 
Malhotra et al. 2001), with the best fit to this sample given by the dotted line (see main text for details). Also shown are 240 2 ; ~ 0 
(U)LIRGs (purple open circles; Diaz-Santos et al. 2013; Farrah et al. 2013), and 12 2 : ~ 0.3 (U)LIRGs (red triangles; Magdis et al. 2014). 
The high- 2 ; comparison samples include 25 1 < 2 ; < 6 star-forming galaxies (cyan squares; D14) and the dashed line indicate the best fit 
to this sample. Also shown are 16 2 ; ~ 4 — 7 quasars (green triangles; lono et al. 2006; Walter et al. 2009; Wagg et al. 2010; Gallerani 
et al. 2012; Wang et al. 2013; Venemans et al. 2012; Garilli et al. 2013; Willott et al. 2013), and 10 normal star-forming 2 ; ~ 5 — 6 galaxies 
(magenta stars) observed by Gapak et al. (2015). The Galaxy is shown for reference with a gray cross (Pineda et al. 2014). 
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S[cii] measurements rather than the IR-color-corrected 
ones. From Fig. 10 we see that over the range in Ssfr 
0.001 — IMq yr“^kpc“^) spanned by these three 
surveys, the H[cii] — ^sfr relation defined by our sim¬ 
ulations is in excellent agreement with the observations. 
Also, the majority of our simulated 1 kpc x 1 kpc regions 
fall within the observed T^sfr and ranges. A small 

fraction (a few procent) of regions in our simulations ex¬ 
hibit an excess in S[cii] for a given Ssfr relative to the 
observed relations, but the majority coincide with the ob¬ 
servations. We note that the observed relations exhibit 
significant scatter 0.2 — 0.3 dex; D14, H15, K15) as 
well as small systematic offsets relative to each other (in 
particular in the case of D14). The S[cii] — Ssfr rela¬ 
tions observed in the five fields in M31 by K15 yield 
best-fit power-law slopes in the range ^ 0.67 — 1.03, 
with an average of 0.77 (shown as dashed-dotted line 
in Fig. 10). Similar power-law fits to the samples of 
D14 and H15 yield slopes of ^ 1.07 and ^ 0.88, re¬ 
spectively (shown as dashed and dotted lines in Fig. 
10). In comparison, a power-law fit to our simulations - 
across the full HsFR-range - results in a slope of ^ 0.60, 
i.e., on the low side of the observed range. However, 
if instead we fit only to simulated regions within the 
^SFR = 0 . 001 — 0.1 Mq yr“^ kpc“^ range, thereby match¬ 
ing the range with the most observations, we find a slope 
of ^ 0.75, i.e., well within the observed range. 

Beyond the above HsFR-range a comparison with the 
observations has to rely on extrapolations of the simple 
power-law fits to the observed il[cii] — ^sfr relations. 
At low SsFR (^ 0.001 Mq yr“^kpc“^), the simulations 
broadly follow the extrapolations of the observed rela¬ 
tions, except for D14 where the systematic offset noted at 
higher Hsfr is compounded at the lower Hsfr values ow¬ 
ing to the relatively steep slope of the fit to the D14 data. 
At these low Ssfr levels a larger (but still minor, overall) 
fraction of the simulated regions display excess [Cll] emis¬ 
sion levels (> 10 x) compared to the survey-based power- 
law fits. This excess [Cll] emission is driven by the PDR 
gas in our simulations (second panel in Fig. 10), which 
dominates the [Cll] emission at these Hsfr levels. Inter¬ 
estingly, in at least two of the five fields in M31 studied 
by K15, a similar [Cll]-excess relative to the power-law 
fit is observed at Hsfr ^ 0.001 Mq yr“^ kpc“^ (see Fig. 
7inK15). K15 argues that the excess is due to a con¬ 
tribution from diffuse, ionized gas (Hll), although they 
cannot discount the possibility that the larger dispersion 
is at least partly due to being close to the sensitivity limit 
of their survey at such low Il[cii]- Our simulations, how¬ 
ever, clearly show no significant contribution from the 
diffuse Hll gas phase at low Ssfr- 

At high SsFR (^ 1 Mq yr“^ kpc“^), the scatter in the 
S[cii] — SsFR relation from our simulations decreases sig¬ 
nificantly, and is seen to broadly match the extrapolated 
locally observed relations. Furthermore, the simulation 
contours are seen to agree with the few galaxies ob¬ 
served to date to have Ssfr ^ 1 Mq yr“^ kpc“^: a few 
sources from the D14 sample (there are two D14 galax¬ 
ies with HsFR ^ lA[ 0 yr“^ kpc“^) and ALESS73.1, 
8i z = 4.76 sub-millimeter selected galaxy, marginally 
resolved in [Cll] and with a disk-averaged Ssfr of ^ 
80 Mq yr“^ kpc“^ (De Breuck et al. 2014; indicated by a 
red cross in Fig. 10). 


The aforementioned studies of local galaxies typi¬ 
cally measure both obscured (e.g. from 24 /im) and un¬ 
obscured (from either FUV or Ha) SFRs in order to esti¬ 
mate the total SsFR and, while intrinsic uncertainties are 
inherent in the empirical 24/im/FUV/Ha ^ SFR cali¬ 
brations, the Esfr values should be directly comparable 
to those from our simulations. Nonetheless, some caution 
is called for when comparing our simulations to resolved 
observations, primarily regarding the determination of 
Esfr on <kpc-scales. In particular, as pointed out by 
D14, the emission from old stellar populations in diffuse 
regions with no ongoing star formation can lead to over¬ 
estimates of the SFR when using the above empirical 
SFR calibrations. Both K15 and H15 account for this by 
subtracting the expected cirrus emission from old stars in 
24 /im (using the method presented in Leroy et al. 2012 ), 
although on scales ^ 1 kpc this correction for 24 /im cir¬ 
rus emission becomes particularly challenging as it was 
calibrated on > 1 kpc scales. K15 further points to the 
problem arising from the possibility of photons leaking 
between neighboring stellar populations, meaning that 
average estimates of the SFR on scales < 50 pc might 
not represent the true underlying SFR. 

6.4. Physical underpinnings of the [Cll] — SFR 
relationship 

The SFRs of galaxies at both low and high redshifts is 
observed to strongly correlate with their molecular gas 
content (typically traced by CO or dust emission) (e.g. 
Kennicutt 1998; Daddi et al. 2010; Genzel et al. 2010; 
Narayanan et al. 2012). This SFR — H 2 dependency is 
incorporated into our SPH simulations (see Section 3), 
and it is therefore natural to ask whether the integrated 
[Cll] — SFR relation examined in Section 6.2 might sim¬ 
ply be a case of galaxies with higher SFRs having larger 
(molecular) gas masses with higher associated [Cll] lu¬ 
minosities. To investigate this, we show in Fig. 11 [Cll] 
luminosity vs. gas mass for the full ISM in our simu¬ 
lated galaxies and for the individual gas phases sepa¬ 
rately. For the molecular gas phase we see a significant 
luminosity—mass scaling (red curve), which would ex¬ 
plain the strong L^cu] ~ SFR relation for this phase and, 
at least in part, the L\cu] — SFR relation for the full 
ISM. The PDR gas and the ionized gas do not exhibit 
similar luminosity—mass scaling relations. Since, on av¬ 
erage, significantly more mass resides in each of these two 
phases than in the molecular phase, this has the effect of 
weakening the correlation between T[cn] and SFR when 
considering the full ISM. 

The abundances of metals in the ISM is expected 
to affect the [Cll] emission. In their study of low- 
metallicity dwarf galaxies, D14 compared total (IR+UV) 
SFRs to SFRs derived using the best-fit L^cu] ~ SFR re¬ 
lation to their entire sample. They found an increase in 
the SFRir 0 uv/SFRcii fraction, corresponding to weaker 
[Cll] emission, toward lower metallicities. In Fig. 12 
we show SFR/I/[cii] as a function of ^[o/h] (= 12 + 
log[0/H]) for our simulated galaxies as well as for the 
dwarf galaxy sample of D14. As D14 points out, ^[o/h] 
does not take into account potential deficits of carbon rel¬ 
ative to the [0/H] abundance, potentially obscuring any 
potential correlation between [Cll] luminosity and actual 
carbon abundance. In Fig. 12 we therefore also plot 
SFR/I/[cii] as a function of ^[c/h] (= 12 + log [C/H]) for 
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Fig. 10.— [Cll] luminosity surface density (S[qjj]) vs. SFR surface density (Ssfr) for the full ISM in our simulated galaxies (top panel), 
and for each of the three ISM phases (bottom three panels), '^[cu] ^^^<1 ^SFR are determined over 1 kpc X 1 kpc regions within all 7 galaxies. 
The filled colored contours indicate the number of such regions with this combination of and as a percentage (at 0.1%, 1%, 

10%, 40% and 70%) of the maximum number of regions. The solid line shown in all panels is the best-fit power-law to the total gas '^[cu] 
vs. Ssfr (see legend). For comparison we show observed vs. Ssfr for individual 20pc regions in M31 (green dots; K15), 1 kpc 

regions in nearby (mostly spiral) galaxies (magenta contours, indicating regions containing 25%, 45%, and 95% of the total number of data 
points; H15), and 1 kpc regions in local dwarf galaxies (cyan crosses; D14). The best-fit power-laws to these three data-sets are shown as 
dashed-dotted, dotted, and dashed lines, respectively. In the case of M31 the power-law is fitted to 50 pc regions (see K15). We also show 
the galaxy-averaged S[cii] and Esfr of ALESS73.1 (red cross), a, z = 4.76 submillimeter-selected galaxy which was marginally resolved in 
[Cll] and in the FIR continuum with ALMA, revealing a source size of R ~ 2 kpc, = 5.15 x 10^ Lq, and SFR ~ 1000 M© yr“^ (De 

Breuck et al. 2014). 
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Fig. 11.— [Cii] luminosity versus gas mass of the total ISM (black 
line), the molecular (red line), the PDR (orange line), and the Hll 
(blue line) gas phase for each of our seven simulated galaxies. The 
[Cll] emission from the molecular gas is seen to increase with the 
amount of molecular gas available, whereas this is not the case for 
gas associated with PDR and Hii regions. 

our simulated galaxies. Arguably, our simulated galax¬ 
ies exhibit a weak trend in SFR/I/[cii] with metallicity 
(in particular for [C/H]). Certainly, there is a very good 
overall agreement with the findings of D14. The lack of a 
strong trend within our simulation sample is not supris¬ 
ing, however, since our galaxies span a limited range in 
metallicity (^[o/h] = 8-03 — 8.96) and, as we saw in Sec¬ 
tion 6.2, form a tight Lcii — SFR relation. 
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Fig. 12.— SFR/L[qjj] as a function of metallicity for our simu¬ 
lated galaxies. The metallicity is parametrized both relative to the 
oxygen (circles) and the carbon (crosses) abundance. The galaxies 
are color coded in the usual manner. For comparison, we show 
SFR/L[cii] vs. ^[o/H] tor the local dwarf galaxy sample of D14, 
where SFRs are from combined UV and 24 p,m measurements. 


To better understand how the metallicity and also the 
ISM pressure affects the [Cll] emission, we define a ‘ [Cll] 
emission efficiency’ (e[cii]) of a given gas phase as the 
[Cll] luminosity of the gas phase divided by its mass, and 
examine how it depends on Z' and Pext /• Specifically, 
we create a grid of (Z', Pext/^e) values, and in each 
grid point the median e[cii] of all clouds in our simulated 
galaxies is calculated. The resulting e[cii] contours as 
a function of Z' and Pext/^B are shown in Fig. 13 for 
the molecular (top), PDR (middle), and Hll (bottom) 


phases. For comparison, we also show the median Z'- 
and Pext/^B-values for each of our galaxies. 
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Fig. 13.— Contours of (median) e[cii] (he. [Cll] luminosity per 
gas mass) as a function of Z' and Pext / for the molecular (top), 
PDR (middle), and Hll (bottom) gas phases in our model galaxies. 
In each panel the contours indicate 0.01%, 1%, 2%, 5%, 10%, 20%, 
50%, 70%, and 90% percentages of the maximum efficiency. For 
reference, the maximum (median) e[cii] for the molecular, PDR, 
and ionized gas phases are: 3.7, 6.6, and O.O2L0 /Mq, respectively. 
Median values of Z' and Pext/^B and for the CMC and ionized 
cloud population in each model galaxies are indicated as colored 
circles (error bars are the Icr dispersion of the distributions). 


The molecular phase in our simulations is seen to radi¬ 
ate most efficiently in [Cll] at relatively high metallicities 
(logZ' ^ 0.5) and cloud external pressures (Pext/^B ^ 
10^cm“^K). The PDR and ionized phases, however. 
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have their maximal e[cii] at \ogZ' > 0 and Pgxt ^ 
10^cm“^K, i.e. at ^ 100x lower pressures than the 
molecular phase, and can maintain a significant [Cll] ef¬ 
ficiency (> 20%) even at Pext/^B ^ 10^ cm“^ K provided 
logZ' ^ 0. This is consistent with our finding in Sec¬ 
tion 6.1, that the molecular phase dominates the [Cll] 
emission in the central, high-pressure regions, while the 
PDR gas dominates further out in the galaxy where the 
ISM pressure is less extreme. In all three panels, G6 and 
G7 lie in regions of high e[cii] 50% for the molecular 
and PDR gas). Thus, the reason why these two galaxies 
have the highest total [Gll] luminosities is in part due to 
the fact that they have the highest molecular gas masses 
(by nearly an order of magnitude, see Fig. 11), and in 
part due to the bulk of their cloud population (GMCs or 
ionized clouds) having sufficiently high metallicities (on 
average > 5x higher than G2) and experiencing external 
pressures (^ 100 x higher on average than the remaining 
galaxies) which drives their [Gll] efficiencies up. 

7. COMPARING WITH OTHER [Cii] SIMULATIONS 

A direct quantitative comparison with other [Gll] emis¬ 
sion simulation studies in the literature (Nagamine et al. 
2006; Vallini et al. 2013, 2015; Popping et al. 2014a; 
Munoz & Furlanetto 2014) is complicated by the fact 
that there is not always overlap in the masses and/or 
redshifts of the galaxies simulated by the aforementioned 
studies and our model galaxies. Furthermore, there are 
fundamental differences in the simulation approach, with 
some adopting semi-analytical models (Popping et al. 
2014a; Munoz & Furlanetto 2014) and others adopting 
SPH simulations (Nagamine et al. 2006; Vallini et al. 
2013, 2015). Also, differences in the numerical resolu¬ 
tion of both types of simulations, and in the specifics of 
the sub-grid physics implemented, can lead to diverging 
results and make comparisons difficult. 

The simulations presented here combine cosmological 
SPH galaxy simulations with a sub-grid treatment of a 
multi-phase ISM that is locally heated by FUV radiation 
and CRs in a manner that depends on the local SFR den¬ 
sity within the galaxies. We have applied our method to 
a high resolution {drnspu — 1.7 x 10^ Mq) cosmologi¬ 
cal SPH simulation of z = 2 star-forming galaxies (i.e., 
baryonic mass resolution dmspYL — 1-7 x 10^ Mq and 
gravitational softening length 0.6h“^kpc, see Section 
3). A novel feature of our simulations is the inclusion of 
molecular, neutral and ionized gas as contributors to the 
[Cll] emission. Another unique feature of our model is 
the inclusion of CRs as a route to produce deep in¬ 
side the GMCs where UV photons cannot penetrate (see 
Section 5). 

Our simulations probably come closest - both in terms 
of methodology and galaxies simulated - to those of 
Nagamine et al. (2006) who employed cosmological simu¬ 
lations with GADGET-2 (the precursor for GADGET- 
3 used here) to predict total [Cll] luminosities from dark 
matter halos associated with z cz 3 LBGs with ^ 
10^^ Mq and SFR > 3OM0yr“^ (see also Nagamine 
et al. 2004). Their simulations employed gas mass reso¬ 
lutions in the range Smspu ^ 3 x 10^“^ Mq and gravita¬ 
tional softening lengths of typically > lh“^kpc. In their 
simulations the ISM consist of a CNM (T ^ 80K and 
n 10cm“^) and a warm neutral medium (T ^ 8000 K 


and n ^ 0.1 cm“^) in pressure-equilibrium, and the as¬ 
sumption is made that the [Cll] emission only originates 
from the former phase. The thermal balance calculation 
of the ISM includes heating by grain photo-electric effect, 
CRs, X-rays, and photo-ionization of C l - all of which are 
assumed to scale on the local SFR surface density. Their 
simulations predict r\j (0.3-1) X 108 L© for their 

brightest LBGs (SFR ^ 30 Mq yr“^) which is slightly be¬ 
low the prediction of our integrated [Cll] — SFR relation 
(^[cii] ^ 6 X 10^ Lq; Fig. 9). 

Employing a semi-analytical model of galaxy forma¬ 
tion, Popping et al. (2014a) made predictions of the in¬ 
tegrated [Cll] luminosity from galaxies at ^ = 2 with 
M* - 10^ - 10^2 M© and SFR - 0.1 - 100 M© yr-^ (see 
also Popping et al. 2014b). The galaxies are assumed 
to have exponential disk gas density profiles, with ran¬ 
domly placed ‘over-densities’ mimicking GMCs (consti¬ 
tuting ^ 1% of the total volume). The gas is embedded 
in a background UV radiation field of 1 Habing, with lo¬ 
cal variations in the radiation field set to scale with the 
local SFR surface density. The C+ abundance is set to 
scale with the abundance of carbon in the cold gas. The 
excitation of [Cll] occurs via collisions with e“, H, and 
H 2 , and the line emission is calculated with a 3D radia¬ 
tive transfer code that takes into account the kinematics 
and optical depth effects of the gas. For galaxies with 
SFRs similar to our model galaxies (^ 5 — 50M©yr“^) 
their simulations predict ensemble-median [Cll] luminosi¬ 
ties of ^ (1 — 6) X 10^ L©.^ This is lower than the [Cll] 
luminosities predicted by our simulations and also some¬ 
what on the low side of the observed 2 : r\j 0 ^[cii] ~ SFR 
relation (Fig. 9). Allowing for the +2(7 deviation from 
the median of the Popping et al. (2014a) models results 
in L[cii] < (2 - 30) X 10^ L© for SFR - 5 - 50 M© yr-^, 
which matches the observed L^qu] ~ SFR relation. 

The remaining [Cll] simulation studies in the litera¬ 
ture focus on ^ ^ 6 galaxies (Vallini et al. 2013, 2015; 
Munoz & Furlanetto 2014). Vallini et al. (2013) uses 
a GADGET-2 cosmological SPH simulation with a mass 
resolution drnspu = 1.32 x 10^ M© and gravitational soft¬ 
ening length ^ 2h~^ kpc. They adopt a two-phased ISM 
model (GNM+WNM), and heating and cooling mecha¬ 
nisms similar to that of Nagamine et al. (2006), in or¬ 
der to predict the [Gll] emission from a 2 : = 6.6 Lyman 
alpha emitter (LAE) with SFR 10M©yr“^. They in¬ 
vestigated two cases of fixed metallicity: Z' = 1 and 
Z' = 0.02. In their simulations, the CNM (T ^ 250 K 
and n ^ 50 cm“^) is found to be responsible for ^ 95% of 
the total [Cll] emission, with the WNM phase (T ^ 500 K 
and n ^ 1 cm“^) contributing the remaining 5%. Vallini 
et al. (2015) presents an update to their 2013 model, in 
which the same ^ = 6.6 SPH simulation as in Vallini 
et al. (2013) is considered but now with the implementa¬ 
tion of a density-dependent prescription for the metallic- 
ity of the gas, and the inclusion of [Cll] contributions 
from PDRs (in addition to their previous two-phased 

® Since Popping et al. (2014a) plots against Ljk and not 

SFR, we have converted our SFRs to Ljk in order to crudely esti¬ 
mate their predicted [Cll] luminosities (see their Fig. II). For the 
SFR —)► Liji conversion we have used (Bell 2003). Not all the star 
formation will be obscured and so the IR luminosities, and thereby 
the [Cll] luminosities given here will be upper limits. 
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CNM+WNM ISM model), and accounting for the ef¬ 
fect of the CMB on the [Cll] emission. As a result of 
these updates, it is found that the [Cll] emission is now 
dominated by the PDRs, with < 10% coming from the 
CNM. This is qualitatively consistent with the results 
from our simulations at 2: = 2 where the PDR gas domi¬ 
nates the total [Cll] emission at least at the low SFR end 
10 Mq yr“^). However, our simulations do not incor¬ 
porate the CNM to the same extent as that of Vallini 
et al. (2015), as the lowest densities found in the neutral 
gas in our simulations is of order ^ 100 cm“^, i.e., above 
typical CNM densities of ^ 20 — 50cm“^. It is therefore 
reassuring that Vallini et al. (2015) find the CNM contri¬ 
bution to the total [Cll] emission to be benign. Assuming 
that HsFR oc I;h 2 and I;h 2 H[cii], Vallini et al. (2015) 
scale the [Cll] luminosity of their fiducial LAE model 
(SFR = 10 Mq yr“^) in order to generate a L^cu] ~ SFR 
relation. 

Munoz & Furlanetto (2014) make analytical predic¬ 
tions of the [Cll] luminosities for a range of galaxy types 
at ^ ^ 6 (e.g., Lyman-alpha emitters, starburst galax¬ 
ies and quasars, spanning a range in SFR from tens of 
MQyr“^ to several thousand) as part of their efforts to 
develop an analytical framework for disk galaxy forma¬ 
tion and evolution at these early epochs. In their study, 
the [Cll] emitting gas is assumed to come from photo¬ 
dissociation regions only. Throughout their models, the 
metallicity is kept fixed at solar. By tuning their mod¬ 
els, i.e., either increasing the star formation efficiency 
at high redshifts or lowering the depletion of metals 
onto dust grains, they arrive at the [Cll] —SFR relation: 

L[cii]/Lo = 5 X 10® (SFR/100M©yr-i)° t While this 
is essentially a linear relation, it struggles to match the 
expected [Cll] luminosities based on observations due to 
the somewhat lower normalization. 

8. CONCLUSION 

We have developed SIGAME to include simulations of 
the [Cll] emission from star-forming galaxies. The code 
employs a multi-phased ISM consisting of UV- and CR- 
heated clouds of molecular and PDR gas as well as diffuse 
regions of ionized gas, and traces the [Cll] emission from 
each of these three phases. SIGAME was applied to SPH 
simulations of seven star-forming galaxies at 2: = 2 with 
stellar masses in the range r\j (0.4 — 6.6) X 10^^ Mq and 
SFRs ^ 5 — 60MQyr“^ in order to make predictions of 
the [Cll] line emission from MS galaxies during the peak 
of the cosmic star formation history. 

A key result of our simulations is that the total [Cll] 
emission budget from our galaxies is dominated by the 
molecular gas phase (^ 70%) in the central regions 
{R ^ Ikpc) where the bulk of the star formation oc¬ 
curs and is most intense, and by PDR regions further 
out (R > 1 — 2 kpc) where the molecular [Cll] emission 
has dropped by at least an order of magnitude compared 
to their central values. The PDR gas phase, while rarely 
able to produce [Cll] emission as intense as the molecular 
gas in the central regions, is nonetheless able to maintain 
significant levels of [Cll] emission from R ^ 2 kpc all the 
way out to ^ 8 kpc from the center. The net effect of this 
is that on global scales the PDR gas can produce between 
8% and 67% of the total [Cll] luminosity with the molec¬ 
ular gas responsible for the remaining emission. We see 


a trend in which galaxies with higher SFRs also have a 
higher fraction of their total [Cll] luminosity coming from 
the molecular phase. Our simulations consistently show 
that the ionized gas contribution to the [Cll] luminos¬ 
ity is negligible (< 3%), despite the fact that this phase 
dominates the ISM mass budget (see Fig. 11). Therefore, 
the ionized gas phase is an inefficient [Cll] line emitter in 
our simulations. 

The integrated [Cll] luminosities of our simulated 
galaxies strongly correlate with their SFRs, and in a 
manner that agrees well (both in terms of slope and 
overall normalization) with the observed L^cu^ — SFR 
relations for normal star-forming galaxies at both low 
and high redshifts. We have also examined the rela¬ 
tionship between the 1 kpc-averaged surface densities of 
I/[cii] and SFR across our simulated galaxies. The re¬ 
sulting — Fsfr relation spans six orders of mag¬ 

nitude in Fsfr (^ 10“^ — 10 Mq yr“^ kpc“^), extend¬ 
ing beyond the observed ranges at both the low and 
high end of the relation. In the FsFR-range where a 
direct comparison with observations can be made (^ 
0.001 — 1 Mq yr“^ kpc“^) we find excellent agreement 
with our simulations. 

Our simulations suggest that the correlation between 
[Cll] and SFR - both the integrated and the resolved ver¬ 
sions - is determined by the combined [Cll]-contribution 
from the molecular and PDR phases (the ionized gas 
makes a negligible contribution), with the former exhibit¬ 
ing the steepest slope and dominating the [Cll] emission 
at the high-SFR-end. We argue that this is due to the 
fact that the [Cll] luminosity scales with the amount 
of molecular gas present in our simulated galaxies. A 
similar luminosity-mass scaling is not seen for the other 
phases. Our work therefore suggest that the observed 
[Cll] — SFR relation is a combination of the line predom¬ 
inantly tracing the molecular gas (i.e., the star forma¬ 
tion ’fuel’) at high SFR levels/surface densities, while 
at low SFRs/surface densities the line is tracing PDR 
gas being exposed to a weaker, interstellar UV-field. As 
a consequence, we hypothesize that galaxies with large 
mid-plane pressures and large molecular gas fractions will 
display a steeper [Cll] — SFR relationship than galaxies 
where a larger fraction of the ISM is atomic/ionized gas. 
In the future we will extend this study to a larger sample 
of model galaxies, in particular with a larger spread in 
SFRs and metallicities. 
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APPENDIX 

A. [Cii] EXCITATION AND EMISSION 


TABLE 2 

Provides a summary of the various heating and cooling mechanisms adopted in the molecular and atomic gas regions of 

OUR CMC MODELS, ALONG WITH THE PHYSICAL QUANTITIES ON WHICH THEY DEPEND. 


Process 

Cooling and Heating Rates in CMC Models 
Parameters 

Reference 

PpE 

Photo-electric heating 

Go, Tk, Re, Rh 

Bakes & Tielens (1994) 

rCR,Hi 

Cosmic ray heating in atomic gas 

CcR, '^Ho Xe 

Draine (2011) 

rcR,H2 

Cosmic ray heating in molecular gas 

CCR, RH2 

Stabler & Palla (2005) 

^H2 

H2 line cooling 

RH2 5 

Papadopoulos et al. (2014) 

Aqi 

Ol line cooling 

Rh, Tk 

Rollig et al. (2006) 

Acii 

[Cll] line cooling in ionized gas 

Go, Tk, rh, ATc, ay 

Coldsmith et al. (2012) 


Note. The most important references are also given, and we also refer to Olsen et al. (2015) for a detailed description. 


For a two-level system such as [Cll] embedded in a radiation field with energy density U at the transition frequency, 
the general rate equation governing the population levels can be written as: 


ni 


BiuU + Cu 


__ gBuiU gKCui 

^ui + BuiU + Cui Aui + BuiU + Cui 


(Al) 


where Aui {= 2.3 x 10 ^s is the spontaneous emission rate, and Bui and Biu are the stimulated emission and 
absorption rate coefficients, respectively (Goldsmith et al. 2012). Here we have invoked detailed balance, i.e. BiuU = 
gBuiU and Ciu = gKCuh where g (= gu/gi) is the ratio of the statistical weights and K = = g-9i.25K/Tk^ 

We can then write the fraction of in the upper level, fu, as: 

no 1 

fu = 


ni ^Uu 




gBulU + gKCui 


^ui + (1 + g)BuiU + Cui + gKCui 
_ gK -h gBuiU/Cui _ 

1 + gK + AuijCui + (1 + g) BuiU /Cui 


(A2) 


In general, U = {1 — /3)U{Tq^) + /3[/(Tbg), where U{T\^g) is the energy density from a background field (e.g. CMB or 
radiation from dust), and P is the escape probability fraction at the [Cll] frequency (see also Goldsmith et al. 2012). 
We shall ignore any background, however, in which case we can write the stimulated downward rate as: 


and the excitation temperature as: 


BuiU = 


-91.25 K/T®^ 


(1 - ^)Aul 


K 


-U 

(A3) 

l3Aui'\ 

Cui ) ’ 

(A4) 


see Goldsmith et al. (2012). For the escape probability we assume a spherical geometry with a radial velocity gradient 
proportional to radius, such that /3 = (1 — exp(—r))/r. In calculating the integral in eq. 18 we split the integral up 
into 100 radial bins of thickness AR (see Section 5), and we approximate the optical depth in each such bin with that 
of a homogeneous static slab of gas of thickness AR (Draine 2011): 


T = 


gu Auic 
gi 4(27r)3/2z/3cr^ 


niAR 


L _ nugi \ 

V ^iguj ’ 


(A5) 


where ay is the gas velocity dispersion, which is calculated according to eq. 5 for the PDR and molecular gas regions 
in the GMGs, and is set to the local velocity dispersion in the SPH simulation in the case of the ionized clouds. We 
use eqs. A2, A3, A4, and A5 to iteratively solve for consistent values of fu and (3 in each AR bin. This is done by first 
assuming optically thin emission (/3 = 1) in order to get an initial estimate of fu (eq. A2), which is subsequently used 
to calculate r and (3 (eq. A5) and from that Tex and BuiU (eqs. A4 and A3), etc. Once consistent values for fu and r 
have been reached, we calculate the total cooling rate according to eq. 19. This cooling rate is used to determine the 
thermal balance at various points within the GMGs as well as the [Gll] emission from the ionized clouds as described 
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in Section 5. We emphasize that our methods assumes that the [Cll] emission from the different AR bins within a 
cloud is radiatively de-coupled, and that the total [Cll] emission from a cloud is therefore the sum of the emission from 
each bin. 

As explained in Section 5 [Cll] is collisionally (de)excited by H2 in the molecular phase, by e~ and Hi in the PDR 
region, and by e~ in the ionized gas. For a single collision partner the collision rates are equal to the density (n) of the 
collision partner times the rate coefficients, i.e., Cui = nRui and Ciu = nRiu- In case of two collision partners we have 
Cui = and Ciu = + ^2^Zu,2- Fig- 14 shows the [Cll] deexcitation rate coefficients that we 

have adopted in our work for collisions with e~ Hi, and H2 as a function of temperature. For collisions with electrons, 
we adopt the following expression for the deexcitation rate coefficient as a function of electron temperature (Te): 

Rui{e-) = 8.7 X 10-^(Te/2000K)-^-^^cm^s-\ (A6) 

which is applicable for temperatures from 100 K to 20,000 K (Goldsmith et al. 2012). At temperatures > 20,000 K we 
set Rui{e~) (Wilson & Bell 2002; Langer & Pineda 2015). For the PDR gas, the electron density is calculated 

using GLOUDY and the electron temperature is set to the kinetic temperature of the gas (calculated according to eq. 
16). For the ionized gas we assume Ue = nnn, i.e., Xe = 1, and the electron temperature gas is set to the SPH gas 
kinetic temperature (see Section 4.2). For the deexcitation rate coefficient for collisions with Hi we use the analytical 
expression provided by Barinovs et al. (2005) and shown as the dotted curve in Fig. 14: 

= 7.938 X 10^^ exp(-91.2/Tk) (l6 + 0.34477^ - 47.7/Tk) cm^ s“^ (A7) 

While Barinovs et al. (2005) cites an application range for the above expression of 15K < < 2000 K, a comparison 

with (Hi)- values found by Keenan et al. (1986) over the temperature range lOK < Tk < 100, 000 K shows that 

eq. A7 provides an excellent match over this larger temperature range (Fig. 14). For collisions with H2 we follow 
Goldsmith et al. (2012) and assume that the collision rate coefficients are approximately half those for collisions with 
H over the relevant temperature range for the molecular gas, i.e. Rui(ii 2 ) = 0.5 x 



log( T [K] ) 

Fig. 14.— [Cll] deexcitation coefficients {Rul) as a function of temperature for collisions with e~ (dashed line), Hi (dotted line), and 
H 2 (solid line and diamonds). The curves are thicker over the temperature ranges where they are formally applicable. The blue, orange, 
and red rectangles indicate the temperature range encountered in the ionized, atomic, and molecular phase in our simulated galaxies, 
respectively (see sections 4.1.2 and 4.2). 

















